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Ultracold atoms in optical lattices undergo a quantum phase transition from a superfluid to a 
Mott insulator as the lattice potential depth is increased. We describe an approximate theory of 
interacting bosons in optical lattices which provides a qualitative description of both superfluid and 
insulator states. The theory is based on a change of variables in which the boson coherent state 
amplitude is replaced by an effective potential which promotes phase coherence between different 
number states on each lattice site. It is illustrated here by applying it to uniform and fully frustrated 
lattice cases, but is simple enough that it can easily be applied to spatially inhomogeneous lattice 
systems. 
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I. INTRODUCTION 

The observation [1] of a cold-atom quantum phase 
transition between superfluid (SF) and Mott insulator 
(MI) states was important on its own merits, and also 
because it suggested future experimental studies of clean 
highly controllable strongly correlated bosonic many- 
body systems. The promise of early experiments has 
been borne out by studies that have demonstrated un- 
precedented experimental control in designing and in- 
vestigating many body systems whose Hamiltonian's are 
known with a level of precision that is uncommon in 
condensed matter physics [2j. Cold atom systems are 
not, however, completely free of the real world complica- 
tions that can confuse the interpretation of experiments. 
The most obvious troublesome complications in simulat- 
ing condensed-matter many-body physics problems us- 
ing cold atoms are that experimental systems are always 
spatially inhomogeneous to some degree, and that they 
are often fairly small. In most cases, the spatial inho- 
mogeneity is an undesirable consequence of an experi- 
mental necessity, for example the harmonic trapping po- 
tential employed in most cold atoms set-ups. In some 
cases, though, it is the central focus of the experiment, 
as in studies of disorder in strongly correlated bosonic 
systems [H, 0] • In this paper we describe an approximate 
theory of strongly interacting bosons in an optical lattice 
that is sophisticated enough to achieve a good qualitative 
description of both Mott insulator and superfluid limits 
and yet simple enough that it can be applied with relative 
ease to finite spatially inhomogeneous bosonic optical lat- 
tice experiments. The theory is a generalization of the 
mean-field theory of the MI-SF phase transition in which 
the potential which induces coherence between different 
number states on a given lattice site is elevated from a 
variational parameter to a fluctuating quantum variable. 
We illustrate the potential of this simple theory by ap- 
plying it to uniform optical lattice systems with constant 
and fully frustrated inter-site hopping parameters. 

The systems in which we arc interested provide an ap- 



proximate realization of the Bose Hubbard Hamiltonian 
(BHH) 1, 
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The BHH provides an accurate description of cold atom 
systems in which the optical lattice potential is strong 
enough that only the lowest Bloch band is significantly 
occupied 0. In Eq. ((T|), a\ is the boson creation oper- 
ator at site i, rii = a\a i is the number operator, f,-j is 
the hopping amplitude between sites i and j, U is the 
on-site interaction energy, fj, is the chemical potential, 
and ei is an energy offset due to the trap, or to other 
intended or unintended local potentials. For a transla- 
tionally invariant system with nearest-neighbor hopping, 
mean field theory produces a phase diagram in fi/U-t/U 
space in which SF states are interrupted at small t/U by 
a series of MI lobes centered on half-odd integer values 
of [i/U , each characterized by a different fixed integer 
value N of the number of atoms per site (Fig. [T]). Since 
t decreases and U increases with optical lattice potential 
strength, t/U can be experimentally varied over a wide 
range. 

Most studies of the BHH have focused on some regime 
of the model's parameter space. For instance, in the 
small t/U limit atom number fluctuations on a given site 
due to hopping can be treated as weak perturbations. 
Even near the SF-MI transition strong interactions still 
suppress number fluctuations significantly, reducing the 
physically relevant Fock subspace to two or three num- 
ber states and justifying mappings which transform the 
BHH into spin models that can be attacked using a large 
arsenal of extensively developed techniques [7|. On the 
other hand, for large values of t/U the interactions be- 
tween cold atoms are weak enough to justify Bogoliubov's 
weakly interacting boson theory. For large numbers of 
atoms per unit cell one can often employ the rotor ap- 
proximation, a, = \/fle iei that is valid when the mean 



2 



3.0 
2.5 
2.0 
,1.5 
1.0 
0.5 



0.0 
0.00 





— 


.... 


^T"^ 



















0.05 



0.10 



0.15 

zt/U 



0.20 



0.25 



0.30 



FIG. 1. Mean field phase diagram of the BHH. The solid 
lobes correspond to the MI phases, characterized by integer 
occupation of atomic sites, z is the coordination number, 
the number of neighbors of any given site. For sufficiently 
large values of zt/U, the system enters the SF phase. The 
dashed curves are contours along which the coherence field 
Berry curvature vanishes (see text). 



occupation number n at each site is so large that its rel- 
ative fluctuations are small. The resulting Hamiltonian 
is a quantum phase model in which the degrees of free- 
dom are the phases of the superfluid at different sites. In 
this limit interactions induce phase fluctuations around 
a mean-field state in which all sites adopt a common 
phase iteSachdev. All of these approaches have disad- 
vantages in describing realistic optical lattice experimen- 
tal systems which may have local superfluidity in one 
part of the system and local insulating behavior in an- 
other, and which typically have a mean boson number 
on each lattice site of order 1 [H, [§] . Our approach has 
goals that are similar to those of other complementary 
approximate theories [Iol - [26j . We seek an approach that 
can adequately describe the physics of the BHH in both 
insulating and superfluid regimes and is simple enough 
that it can be applied in the presence of inhomogeneities. 
One advantage of our approach is that we treat interac- 
tions exactly, a feature most useful in describing strongly 
correlated systems. 

Our paper is organized as follows. In Sec. II, we de- 
scribe the flexible formalism that is the subject of this 
paper. In practical applications it leads to a quadratic ac- 
tion for elementary excitations of uniform or non- uniform 
interacting bosons. In Section III, we report on illustra- 
tive applications first to the case of a uniform BHH with 
constant intersite tunneling amplitudes and then, as an 
example of a non-uniform system, to the case of a uniform 
BHH with fully frustrated intersite tunneling amplitudes. 
In Sec. IV we discuss the limitations of our theory before 
concluding with a brief summary. 



II. FORMALISM 

Our approach is based on single-site interacting boson 
wavefunctions (^(E)) which depend on a complex param- 
eter E and are defined as Fock-space normalized ground 
states of the following single-site Hamiltonian, 



/i(E) = -^ n i n ~ 1) — l m ~ — 



(2) 



Notice that the potential £ induces coherence between 
single site states with different boson occupation num- 
bers. The mean field theory of the BHH SF-MI phase 
transition [27j can be derived by considering variational 
wavefunctions of the following form, 



(3) 



These mean-field wavefunctions do not allow for corre- 
lated inter-site fluctuations. The mean-field ground state 
is determined by minimizing 



£(£) 



(y(Z)\H BHH \*(E)) 
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(4) 



with respect to the variational parameter E. In the SF 
state Emf ^ 0. 

Our approach is to elevate £ from a variational param- 
eter to a quantum variable with correlated spatial fluctu- 
ations by allowing it to depend on site and on imaginary 
time (E — > Ej(r)) and then to construct an action S 
which depends on these fluctuations [HI, [29| , 



S= dr 
Jo 



^ i (V;(E i (r))|^(E i (r))) i + J B[E] 



Here at each instant of imaginary time 

= (*[S]|ff Bgg |tt[S]) 
1 J (tf[£]|tf[£]> 

and the correlated product state is given by 
|*[E])=JJ^(E < )) < . 



(5) 
(6) 

(7) 



In practice the action can be evaluated analytically only if 
the coherence fields are expanded to leading order around 
their mean-field values. The main advantage of this ap- 
proach, as stressed above, is its convenience in practical 
calculations, especially for non-uniform systems. Before 
we elaborate on this point, we examine some formal prop- 
erties of our single site interacting boson wavefunctions 
|tf(E)). 



A. Formal Properties of the Wavefunction 

One method of characterizing the Fock space wave- 
functions |-0(E)) is to consider their expansion in terms 
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FIG. 2. The overlap of our interacting wavefunctions with Fock number states, \(n\tp(T,MF))\ 2 ■ Deep inside the SF phase (large 
|Smf|), more number states come into play, the result of strong atom number fluctuations. At /x = OAU, only one Fock state, 
viz. \n — 1), is dominant close to the phase boundary (small |Emf|), whereas there are two, |n = 1) and \n — 2), at fi = U, in 
accordance with Fig. [1] 



of number eigenstates, 



=X> n (|£|)e X p(m<fe)|n) 

71 = 



(8) 



where we have noted that the magnitude of the expansion 
coefficients depends only on and defined 0s as the 



phase of E. In Fig. [2] we plot 



\(n\mW vs. |E| 



for a variety of n values for both fj,/U = 0.4, which falls 
inside the N = 1 MI lobe and for fi/U = 1.0 at the 
boundary of the N = 1 and N = 2 MI lobes. This plot 
illustrates why spin model approximations to the BHH 
arc justified close to the transition (for small values of 
|S|), since a small number of number states dominate. 
As one goes deeper into the SF region, the potential £ 
induces coherence amongst more number states and spin- 
model approximations will fail. 

The mean field state is characterized by a time inde- 
pendent field Em_f that minimizes the energy, 

dE 

= 0. (9) 



MF 



(The first term of (0 does not contribute to the action if 
£ is time- independent.) Quantum fluctuations are incor- 
porated by employing a Gaussian-fluctuations approxi- 
mation, Sj = Yimf + Since the time dependence 
comes from the fields only, we expand the time deriva- 
tive in the Berry phase term as 



+ 



d 



(10) 



Substituting this in the action enables us to rewrite the 
first term of (|5|), the Berry phase term, as C^Emf)-^-^, 
where the gauge invariant Berry curvature of site i Ci, 
evaluated at £mf, is given by [30] 
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The Berry phase contribution to the action specifies the 
quantization condition of our fluctuating variables and 
plays an essential role in determining elementary excita- 
tion energies. 

The energy functional, Eq. can also be expanded 
around its mean-field value, 



E[H] = E_ 
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(12) 



where i and j stand for lattice sites and h.c. for complex 
conjugate. Combining the Berry phase term with the 
second order contribution to the energy functional E^ , 
we construct a quadratic action from which we can use 
to calculate the elementary excitations. S = Smf + 

S ( V\z i ,z i ] = / <ItY \C i (E M F)ziZ i + E^\n . (13) 



Zi,Zi], where 



^ [CiCEMF^iZi + EWp 

i 

E^ 2 * [£] contains all the second order terms of ([I 
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FIG. 3. A log-log plot of the Berry curvature as a function 
of the coherence field |Ea/f|. Note the power law decay rule 
d(T,MF) ^ \Y,mf\" a ^ 3 at large |£a/f| values predicted by Eq. 
\VA). This plot is calculated from Eq. (fTTj) . The chemical 
potentials /i are measured in units of U. 



B. Single Site States at Large n. 

Further insight into the properties of the single site 
Hamiltonian, Eq. can be obtained by examining the 
quantum phase model, 

h(E) = —n(n -1)- fin- 2s/n\T,\ cos(<9 - E ), (14) 

that is derived from Eq. © by letting E |E|e ieE , and 
employing the rotor approximation a — > y/ne 18 , which 
is valid when number density fluctuations are small. In 
the rest of this section we write E for |E|. Assuming 
that quantum fluctuations are small allows us to deter- 
mine the average atom number n a ~ (E/[/) 2 / 3 and phase 
9 Q = 9^ by minimizing /i(E) with respect to n and 9. In 
addition, we can expand Eq. (ITTl) to second order about 
the extrema n Q and 9s to arrive at the quadratic Hamil- 
tonian, 



3U 

/i(E) — const H — ~£~( n ~ n °) 2 + 



S 4/3 



(15) 



Phase and atom number are conjugate variables and we 
therefore recognize (|15[) as a quantum harmonic oscillator 
Hamiltonian. Using this analogy we find the energy level 
spacing u = \fiU 1 l 5 Y 1 2 / 3 , mass m = 2/3J7, typical den- 
sity fluctuation Sn = {{n-no) 2 ) 1 / 2 ~ fE/t/) 1 / 3 and typi- 
cal phase fluctuation 89 = ((tf-fls) 2 ) 1 / 2 - (C//E) 1 / 3 . We 
see that density fluctuations are suppressed and phase 
fluctuations enhanced as Emf approaches at the MI 
transition boundary. 



The harmonic oscillator analog, based on the identifi- 
cations p ^ Sn and q f-> 59, can also be used to derive 
an expression for the on-site Berry curvature that is valid 
at large E. The Berry curvature is determined by the 
dependence of the single-site wavefunction on the magni- 
tude and phase of E. We therefore consider the influence 
of perturbations on the eigenstates of Eq. (fT5")) . 



h = — muj q 

2m 2 y 
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X p p - X q q = h - X p p - X q q. (16) 



To first order in A p and X q the ground state of h! is given 
by 



\X p p + X q q\(j> ) 



\<Pn 



(17) 



where e„ and \4> n ) are the eigenvalues and eigenfunctions 
of h. It follows that the harmonic oscillator Berry curva- 
ture is given by (cf. Eq. (fTTj) ) 



C = Im 



d<j> 
dX 2 



<9Ai 



<9Ai 



dX 2 



(18) 



Im 



(<l>o\p\<l>n)(<l>n\q\<l>o) ~ (0ojg|^nK0nWo) 
(^n ^o) 



The matrix elements in this expression can be evaluated 
using [h,p] = imuj 2 q and [h, q] = —ip/m to find 



C(E) = 



1 



(19) 



and hence for our single-site states C ~ E~ 4 / 3 at large 
E. This result is confirmed by the plot (Fig. [3|) of Berry 
curvature values numerically obtained from Eq. (Ill[) . 



III. APPLICATION TO THE BOSE HUBBARD 
HAMILTONIAN 

We now test the theory's practical utility, first by ap- 
plying it to the BHH for a uniform optical lattice, and 
secondly by applying it to the fully frustrated BHH which 
has four atoms per unit cell and therefore introduces in- 
homogeneity. Separating Ei into its mean-field and fluc- 
tuation contributions, the single site Hamiltonian hi be- 
comes 



niijii — 1) * 
hi = [ini - a\T,MF - a-L 



mf - ajzi - a t Zi. 

(20) 



The fluctuations are treated perturbatively. We write 



\^o)i + \lfe)iZ% + \lp£,)iZi, (21) 



where 
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(a)Dispersion inside the MI region. (b)Dispersion at the phase boundary. (c)Dispersion inside the SF phase. 



FIG. 4. Elementary excitation energy ui(k) vs k for a two dimensional uniform BHH along the line, k = k x = k y . Note that 
the spectrum is gapped in the insulating phase, and gapless at and beyond the phase boundary. Also note that for small k, 
the dispersion is linear in the SF region, in accordance with Goldstone's theorem. In this approximate theory the dispersion is 
quadratic at the phase boundary. 



Here \ip n )i and E n are the eigenstates and the energy 
levels of the unperturbed on-site Hamiltonian. With Eq. 
(|2"21 , the Berry curvature can be calculated from Eq. (fTTj) 
quite easily. The energy of the system, given by Eq. ^ , 
can also be expanded (to quadratic order) in terms of the 
fluctuating fields (cf. Eq. fTJ)), 



E = E ( - o) +J2 E ^z l +Y^ (E-z^ZiZj + E^.ZiZjj+h-c-, 

i ij 

(23) 



where we naturally identify 



Ey 



dE 



1 d 2 E 



MF 
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(24) 



MF 



The indices i and j refer to either same site or neighbor- 
ing sites. For a uniform lattice, both the Berry curvature 
and the energy derivatives are independent of site in- 
dices. The mean field state is determined by setting the 
first derivative terms to zero. The MI phase boundary is 
defined by the largest value of t/U for a given [ijU for 
which the energy is miminized by |S| = on all sites. For 
the uniform BHH this procedure reproduces the familiar 
phase diagram plotted in Fig. [1] 



A. Elementary Excitations in a Uniform Lattice 

To determine the elementary excitations of the BHH 
14], we turn to the second order action, Eq. (|13|) . where 
we use Eq. (|23[) for the energy functional. For a uniform 



lattice, we Fourier-transform the fluctuations Zi, 

where /? is the inverse temperature, uj n are the Matsub- 
ara frequencies and N here is the number of sites in the 
lattice. The resulting action is 

S {2) [zz] = — y^\-iu n C{^MF)Zk n Zkn + Z kn A k Z kn 
" k,n 



+ z knBkZ-k-n + z knBkZ-k-n] 



(26) 
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^ k,n>0 



\Zkn Z—k — r 



M(k,n) 



Z-k-T 



where, 
M(k,n) 



-iw n C(T,MF) + Ak Bk 

B k icj n C(T, M F) + A k 



A k = 
B k = 
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d 2 E 



mf d'Eid'Sj 
d 2 E 



mf dY,idY,j 



(cos(k x a) + cos(k y a)), 

MF 

(cos(k x a) + cos(kya)), 
(27) 



MF 



and a is the lattice constant. Here, the indices i and 
j denote neighboring sites. By setting the determinant 
of the matrix M(k, n) in Eq. (|26l) to zero, we obtain an 
expression for the excitation spectrum, 



w(jfe) 



C(Smf) 



(28) 
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FIG. 5. Fully Frustrated Lattice: The thick bonds are the 
links with hopping amplitude ty = —t that frustrate inter- 
site coherence. The dashed box contains a unit cell of the 
frustrated lattice. 



In accordance with Goldstone's theorem, Eq. [28] yields a 
gapless Goldstone mode in the SF phase with linear dis- 
persion at long wavelengths (Fig. 0J;). As the MI phase 
boundary is approached and crossed, excitations become 
more localized and the mode dispersion weakens (Fig. 
Sb) and (Fig. |4ji). At certain points in the phase dia- 
gram the Berry curvature C(T,mf) vanishes (Fig. [TJ and 
our theory of the excitation spectrum becomes unreliable. 
We return to this point in the discussion section. 



B. Fully Frustrated Lattices 

As mentioned above, the theory outlined in Sec. II 
is designed with inhomogeneous systems in mind, and 
is general enough to be applied to systems where the 
translational symmetry is reduced, or altogether lost. We 
demonstrate this with the case of a two-dimensional uni- 
form lattice in which the tunneling amplitude sign alter- 
nates in one direction (See Fig. [5]). This hopping model 
corresponds to a half- flux quantum per square lattice pla- 
quette and is referred to as a fully frustrated BHH [31| . 
For this case we allow translational symmetry to be bro- 
ken by doing the mean-field minimization for a lattice 
with four sites per unit cell. The quadratic theory is sim- 
ilarly modified, with the M matrix in Eq. (|27j) enlarging 
to an 8 by 8 matrix. In general, for the case of an entirely 
inhomogenous lattice with N sites, the mean field state 
is determined by N linearly coupled equations, while the 
elementary excitations are obtained by performing a Bo- 
goliubov diagonalization on the 2N by 2N matrix result- 



ing from the second order action (fl~3|) while preserving the 
bosonic commutation relations [32J, |33j (This procedure 
is discussed in the Appendix.) 

A detailed mean field study of the quantum phase tran- 
sitions in the fully frustrated lattice has been carried out 
elsewhere [H]]. Here we apply Eq. © on each site of 
a plaquette to find the mean field values of the fields, 
EajEbjEc an d Ed. The elementary excitations of the 
system follow from the second order action, Eq. (fT3|) . 
We take advantage of the (reduced) translational sym- 
metry by Fourier decomposing the fluctuation fields at 
each site, 



z s i(t,t) 



1 



I3VN 



k 



z s k{r)e 



ikr. 



(29) 



where I refers to plaquettes and s to the four sites (A, 
B, C and D) within a given plaquette. The Berry phase 
term, calculated at the mean field state defined by the 
above four fields, becomes 



y^^C s z sk (T)d T z sk (T), 



(30) 



while the energy term contains all the second order devia- 
tions of the energy about this mean field state. Replacing 
z sk (t) VC~s z sk i T ) reduces the action into the standard 
form 



\ dr V (z sk {T)d T z sk {T) + E [2) [z sk {T),z sk (T) 
Jo ,. v 



. (31) 

Having elevated the fluctuations into quantum variables 
(obeying bosonic commutation relations), we are now in 
a position to employ the transformation z = Bv, where 

Zsk{r) 
Zs-k{r) 

tions, and perform Bogoluibov diagonalization |32j, |35 
to transform the action into 



that preserves the commutation rela- 



S(2) = / dT Yl (VkdrVk + U)(k)VkVk) , (32) 
Jo k 

from which we can easily identify w(k) as the excitation 
spectrum. A few sample plots of the dispersion are shown 
in Fig. [6] 



IV. DISCUSSION 

As demonstrated above, the Berry phase is the critical 
ingredient in formulating the quantum theory and calcu- 
lating the elementary excitations from it. As such, the 
excitation spectra one obtains become unreliable if and 
when the Berry curvature vanishes. The dashed curves in 
Fig. [T]trace the contours in the phase diagram where the 
Berry curvature, Eq. (jlip , of the uniform BHH becomes 
zero. One point of view towards restoring the quantum 
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FIG. 6. Elementary excitations of the fully frustrated lattice plotted along k = k x = k y . The lowest dispersion is gapped inside 
the insulator, and gapless beyond the phase boundary. Again, obeying Goldstone's theorem, the lowest dispersion in the SF 
phase is gapless, and linear at small k. 



theory is to retain terms that are second order in time 
derivative in deriving the action, Eq. (J5j> [H, HI] . Here, 
we content ourselves to exploring the consequences of a 
vanishing Berry curvature to our theory. To have a bet- 
ter understanding of these special points, we focus on the 
MI phases, where the Berry curvature is explicitly given 
by 



1 



1L 



1 



UJ 



(33) 



For n — 1, for instance, C vanishes at fi/U = y/2 — 1, 
corresponding to a set of points in the phase diagram 
with particle-hole symmetry. In other words, in our for- 
malism, there is no distinction between the definitions for 
the annihilation and creation operators that emerge from 
quantizing the fluctuations, rendering their commutator 
zero. This is consistent with the results of Ref. [5j], where 
quantum phase transitions at the tips of the MI phases, 
the multicritical points with particle-hole symmetry, be- 
long to the universality class of the XY model, where the 
time derivative is second order, in contrast to anywhere 
else across the phase boundary where it is first order, and 
there is an absence of particle- hole symmetry [j| [37J ■ 

In summary, the relative ease with which the above 
theory has determined the basic properties of the fully 
frustrated optical lattice higlights its main focus: inho- 
mogenous systems, such as optical lattices in symmetry 
breaking harmonic traps, or experimental set-ups with 
controlled disorder, which have generated lots of inter- 
est lately (38[. The advantage of our formalism is its 
broad applicability and accuracy despite its simplicity, 
which lends to its convenience. It may also be interest- 
ing to modify the theory developed here for other lattice 



Hamiltonians, such as those involving next neighbor in- 
teractions, where interesting phases such as charge den- 
sity waves and supersolids are predicted jSij]. 
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APPENDIX 

In this section, we discuss some of the diagonalization 
procedure we employed above to arrive at Eq. (|32p . The 
need to perform Bogoliubov diagonalization of many de- 
grees of freedom arises in most considerations of many 
body problems, especially inhomogenous ones, and there 
are now a number of sources in the literature one can 
consult for more details [H, HH . Here we give a more 
succint summary. Let 



z\ 

Zl 
Z N 

z~i 

Z~N 



where N is the total number of sites in the optical lattice. 
As mentioned above, these fluctuations are now quantum 
variables, and thus obey bosonic commutation relations, 
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\z, zH = I, where X = 



with 1 and being 



1 
-1 

the N x N identity and zero matrices, respectively. 

We introduce the Hamiltonian T-L such that second or- 
der energy terms of the action, such as the one in Eq. (f3"l"|) , 
can be written as E^[zi, zi\ = z^T-Lz. Our goal is to per- 
form a canonical transformation z — Bv that diagonalizes 
the Hamiltonian % while preserving the bosonic commu- 
tation relations, 



B^UB = V 



(34) 



where T> = {ei, £2, £2n} are the eigenmodes of interest. 
Using Eq. (l34l) alongside the identities B^IB = I (which 



follows from the commutation relations) and 1 = 1 1 , 
we obtain the expression 



HB = VB, 



(35) 



where 



V = X 1 T> — {ei, £2, . . . , £at, — ejv+i, — £n+2, • ■ • , —^2n}- 

If we consider each vector bi comprising the matrix B, 
we see that Eq. (|3"5|) is essentially an eigenvalue prob- 
lem, where the eigenvalues are particle-hole pairs (e^ = 
— e i+7v), and £j, where 1 < i < N are the fundamental 
excitation modes of the theory such as those plotted in 
Fig. 1 
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